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L.  Introduction 


Consider  the  linear  programming  problem 


T T 

minimize  c x,  + • • • + cR  xR 


subject  to  A..^  x^  = b^ 


A . x,  + • • • + A x =b  s « 2,  K 

si  1 ss  s s 


where 


x is  a vector  of  nt  variables 


matrix 


c.  and  b are  vectors  of  dimension  n and  m respectively, 
t s c s 


n - I V 

t-1 


Such  linear  programs  have  a block  triangular  structure,  as  can 


be  readily  observed  from  the  m x n detached  coefficient  matrix 


was  outlined.  The  central  idea  was  to  use  a square  block  triangular 
'artificial  basis',  that  is,  one  with  square  blocks  on  the  diagonal, 
together  with  a factor  to  correct  for  the  off-square  diagonal  blocks 
of  the  true  basis.  In  essence  this  correction  factor  would  differ 
from  the  m * m identity  matrix  in  as  many  columns  as  the  true  basis 
was  off-square  on  the  diagonal.  The  implication  of  the  persistence 
observation  was  that  these  columns  would  be  small  in  number  and  so 
require  little  work  in  being  maintained  from  one  iteration  to  the  next. 
Further,  to  solve  equations  with  the  square  block  triangular  factor, 
only  the  diagonal  blocks  need  by  factorized.  This  could  yield  a 
substantial  savings  in  storage  requirements,  as  well  as  an  increase 
in  speed. 

In  this  paper  we  shall  extend  these  ideas  by  examining  the 
persistence  property  in  detail,  by  characterizing  factorizations  in 
which  the  correction  factor  has  the  minimal  number  of  structural 
columns,  and  by  showing  that  only  a small  submatrix  of  the  correction 
factor  need  be  maintained  to  execute  the  simplex  algorithm.  Then  we 
shall  develop  in  detail  the  methods  of  updating  the  factorization. 

In  the  remaining  sections  we  shall  discuss  the  implementational  details, 
perform  an  analysis  of  the  work  involved,  and  present  some  computational 
experience. 

In  [14]  Kalllo  and  Porteus  also  employ  a square  block  triangular 
factorization  along  the  lines  suggested  in  [5].  However  their  factori- 
zation and  methods  of  updating  differ  substantially  from  ours.  The 
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factorization  we  consider  always  maintains  the  square  block  triangular 
factor  as  a submatrix  of  A;  in  [14],  this  factor  contains  some  trans- 
formed columns,  making  it  denser  than  ours  and  more  difficult  to  store, 
but  yielding  a sparser  correction  factor.  No  computational  experience 
is  reported  in  [14]. 

Other  methods  devised  specifically  for  this  class  of  problems 
have  considered  mainly  staircase  structures,  that  is,  where  the  matrices 
A = 0 for  t < s - 1.  See  for  example  the  nested  decomposition 
algorithm  of  Ho  and  Manne  [12]  and  the  block  factorization  techniques 
of  Loute  [16]  and  Wollmer  [29].  McBride  [17]  uses  a dynamic  factoriza- 
tion of  the  basis  after  inducing  a block  triangular  structure  with  the 
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use  of  Hellerman  and  Rarick’s  P algorithm  [11].  He  makes  no  prior 
structural  assumptions  and  allows  the  block  triangular  partition  to 
vary  from  one  iteration  to  the  next. 

Bisschop  and  Meeraus  [3]  have  proposed  a method  for  very  general 
L.P.'s  that  is  not  based  on  any  triangularization  or  block  triangulari- 
zation  techniques,  but  uses  an  idea  that  to  some  extent  underlies  the 
method  presented  here  as  well  as  that  of  [17]:  If  we  wish  to  solve 
the  nonsingular  system  Dy  “ d,  but  have  at  our  disposal  a nonsingular 
factorized  matrix  E,  where  E differs  from  D in  (say)  k columns, 
then  there  exists  a suitable  k x k nonsingular  matrix  P such  that 
the  solution  to  Dy  - d can  be  obtained  by  solving  systems  involving 
only  E and  P.  The  method  presented  in  [3]  uses  this  idea  to  save 
substantially  on  storage  costs  over  the  conventional  LU  methods,  by 
confining  the  growth  of  nonzeros  to  P. 


,We  shall  present  empirical  results  that  point  to  the  factoriza- 


tion of  this  paper  being  superior  to  both  the  above  methods  (i.e 
[3],  [17]). 


■MU 


2.  Persistence  in  Staircase  Models 


Time-staged  staircase  linear  programs  are  often  discrete  versions 
of  an  underlying  continuous  process  that  evolves  as  the  solution  to  a 
differential  equation.  One  such  continuous  model  is  the  following 
linear  optimal  control  problem  with  mixed  constraints  on  the  state  and 
control  variables: 


1 m cp  it 

minimize  f {c  x(t)  + d u(t)}dt  + q x(T) 


subject  to  x = Ax  + Bu  + a 


0 = Cx  + Du  + b 


H equations 


x ( t ) > 0 u(t)  >0  t G [0,T]  x(0)  * p 


It  has  been  shown  in  [22]  that  if  (x(*)»  u(-))  is  a nondegenerate 
extreme  point  solution  of  the  above  (convex)  constraint  set,  and  if 
u(-)  is  piecewise  continuous,  then  the  interval  [ 0 , T]  can  be  parti- 
tioned into  a collection  of  open  intervals  (possibly  an  infinite  number) 


such  that  on  each  interval,  precisely  l components  of  (x,u)  are 
strictly  positive,  and  the  remaining  ones  are  identically  zero. 

In  such  cases  we  would  hope  that  as  we  make  finer  and  finer  the 


discrete  approximation  to  the  continuous-time  problem,  the  basic  solu- 
tions will  tend  to  behave  more  and  more  in  this  way,  i.e.,  have  the 


same  activities  basic  over  several  consecutive  time  periods,  and  more- 
over have  precisely  SL  of  them  do  so.  That  this  is  often  the  case 
in  practice  was  observed  by  Dantzig  in  [5]. 

In  the  following  we  shall  investigate  this  persistence  property 
by  looking  at  the  structural  properties  of  bases  in  general  staircase 
models  where  each  step  has  the  same  size, say  Z x r.  Such  models  have 
constraints  of  the  form 


A1  X1  = bl 


Bt-1  Xt-1  + At  Xt  = bt ’ 


> 0 


2,  ...»  K 


where  the  coefficient  matrices  and  B^  are  all  l x r.  To  any 

basis  B will  correspond  matrices  and  that  are  submatrices 


of  and  B^  respectively,  so  that 


*1  P2 


Q2  p3 


qk-i  pk 
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Pj.  and  will  have  the  same  number  of  columns  (which  can  vary  from 

one  t to  the  next);  for  all  t,  Pt  and  Q will  have  £ rows. 

For  t = 1,  K let  kt  be  such  that  Pt  and  Qt  are  £ x (£  + k^)  . 

Clearly  k may  be  of  any  sign  and  is  a measure  of  the  "off-squareness" 
of  period  t in  this  basis.  Note  that  since  B is  square. 


The  following  theorem  gives  bounds  on  how  far  off-square  the 
basis  can  be  along  the  diagonal . 

2.1.  Theorem:  For  all  t and  s such  that  l<_t<^t  + s<^K, 

< l k < * 

j-o  t+j 

Proof:  Consider  the  portion  of  B between  t and  t + s: 
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I 


I 

I 

I 

I 


I 

I 


Since  B is  nonsingular,  the  columns  between  the  two  vertical 
broken  lines  are  linearly  independent.  Therefore,  since  these  columns 
have  only  zeros  outside  the  two  horizontal  solid  lines,  the  number  of 
these  columns  may  not  exceed  the  number  of  rows  between  the  two  solid 
lines,  i.e.. 


I U + k ) < (s  + 2)8. 
j-0  t+J 


or 


l 


j-0 
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The  nonsingularity  of  B also  implies  that  the  rows  between 
the  two  horizontal  broken  lines  are  linearly  independent.  A similar 
argument  yields 


This  is  of  interest  because  it  gives  the  same  bound,  l,  on  the 
off-square  count  for  any  one  period  as  it  does  on  the  count  for  any 
grouping  of  consecutive  periods.  Therefore,  in  a staircase  model 
where  the  number  of  time  periods,  K,  is  very  large  relative  to  £,  we 
can  choose  to  group  together,  say,  every  s periods  to  obtain  a coarser 
partition  of  the  matrix  in  such  a way  that  the  ratio  of  the  off-square 
count  to  the  number  of  periods  in  any  partition,  £./ s , is  small.  Indeed, 
as  K becomes  arbitrarily  large  and  the  number  of  partitions  remains 
constant,  the  off-square  ratio  for  any  partition  becomes  arbitrarily 
small . 

In  many  economic  applications  we  can  attach  a stronger  although 
more  qualitative  significance  to  the  result  of  Theorem  2.1.  Activities 
like  the  level  of  coal  production,  or  the  level  of  an  inventory,  are 
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usually  positive  over  intervals  of  time,  rather  than  at  points  spread 
haphazardly.  In  continuous- time  this  corresponds  to  a statement  about 
the  piecewise  continuity  of  optimal  solutions.  In  discrete-time,  we 
infer  that  such  activities  will  remain  basic  over  a whole  time  interval 
irrespective  of  how  many  time  steps  constitute  the  interval.  The 
implication  of  this  for  the  {kt}  is  that  they  will  be  constant  over 
intervals  of  time  no  matter  how  refined  the  grid  size.  Setting 
kt+j  = k for  j = 0,  1,  ...,  s in  the  theorem  and  noting  that  k is 
an  integer,  we  obtain  k ■ 0 for  grid  sizes  refined  enough,  that 
s > Jl . Thus  in  solutions  where  the  activities  persist  in  the  basis 
over  intervals  of  time,  and  each  such  interval  contains  at  least  2.  + 1 
time  steps,  the  number  of  activities  in  each  time  step  is  precisely 
2.  Note  though,  that  this  will  not  hold  true  at  the  end  points  of  the 
time  intervals. 

In  time-staged  models  where  the  matrix  structure  is  block  tri- 
angular but  not  staircase,  the  above  discussion  and  Theorem  1.1  do  not 
apply.  However,  in  cases  where  the  non-staircase  part  of  the  block 
triangular  structure  contains  just  a sprinkling  of  nonzeros,  and  simi- 
lar type  activities  can  appear  from  one  period  to  the  next,  it  seems 
reasonable  to  expect  a similar  behavior. 
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• The  Square  Block  Triangular  Factorization  of  the  Bas i s 


i 


Let  B denote  the  in  x m basis  matrix.  Since  B is  a sub- 
matrix of  A,  it  too  will  have  a block  triangular  structure 


B - 


J11 


®21  B22 


BK1  BK2  BKK 


Note  that  the  diagonal  blocks  B^^  need  not  be  square  in  general,  and 
that  it  is  possible  that  there  be  no  contribution  to  B from  one  of 
the  'periods',  i.e.,  for  some  i,  the  submatrix 


"kt 


may  be  vacuous.  Note  also  that  since  B is  nonsingular,  it  is  neces- 
sary that  the  leading  submatrices 
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have  full  row  rank  and  hence  have  at  least  as  many  columns  as  rows. 
Likewise  the  submatrices 


BK-1 ,K-1 

bk,k-i  bk,k 


i 


all  have  full  column  rank  and  hence  have  at  least  as  many  rows  as 
columns.  However,  very  little  else  may  be  said  about  B in  general. 

We  shall  factorize  B into  the  product  of  two  m * m nonsingu- 
lar matrices'*- 


where 


B = B F 


(1) 


^Dantzig  [5]  works  with  the  (mathematically)  equivalent  factorization 
B ■»  BE  where  E is  F-*. 


I 
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with  the  diagonal  blocks  Bt(.  square  and  nonsingular.  B will  be 
called  the  artificial  basis. 


To  understand  the  structure  of  F,  consider  the  case  when  the 
i*1*1  column  of  B,  b^,  is  equal  to  the  column  of  B,  b ^ . Letting 

e^  denote  the  jth  column  of  the  m * m identity  matrix,  it  is 


clear  that 


bi  = B ej 


so  that  the  i column  of  F must  be  e ^ . Thus  if  B and  B have 
k columns  in  common  F will  contain  k unit  columns  and  m - k 
other  columns.  By  suitably  permuting  the  rows  and  columns  of  F we 
can  obtain 


F = P 


H I 


where  P and  Q are  permution  matrices,  G is  (m  - k)  x (m  - k) 


and  nonsingular,  and  I is  the  identity  matrix  of  crder  k. 


4.  Minimal  Factorizations 


I 


I 


H 


As  we  shall  emphasize  later  the  usefulness  of  this  factoriza- 
tion approach  will  depend  critically  on  the  number  of  columns  that  B 
and  B have  in  common,  or,  more  directly,  the  size  of  the  matrix  G 
in  (2).  It  will  be  important  to  know  how  to  obtain  a G that  is  as 
small  as  possible. 

4.1.  Definition:  Given  a basis  B,  a square  block  triangular  factori- 
zation of  B,  B = BF,  will  be  said  to  be  minimal  if  for  any  other  such 
factorization,  B = B'F',  B has  at  least  as  many  columns  in  common 
with  B as  does  B'. 

4.2.  Notation:  For  two  matrices  C and  D with  the  same  number 
of  rows,  let  C H D denote  the  submatrix  of  columns  that  are  in 
both  C and  D. 

Thus  a factorization  BF  is  minimal  if  (number  of  columns  of 
B n B)  > (number  of  columns  of  B H B')  for  all  admissable  B' . 

We  have  the  following  useful  characterization  of  minimal 
factorizations. 

4.3.  Theorem:  Let  B be  a given  basis.  Then  a square  block  triangu- 
lar factorization  of  B,  B ■ BF,  is  minimal  if  and  only  if  for 

t ■ 1,  ...»  K,  Btt  H §tt  is  a basis  for  the  space  spanned  by  the 
columns  of  Btt« 
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such  that  the  augmented  matrix 


“ ' [Btt  n i »t> 


has  full  column  rank. 

Since  Btt  is  nonsingular,  bt  has  a representation 


B y = b . 
tt  y t 


Since  W has  full  column  rank  there  is  a column  b of  B, 


0 

b^ 


(3) 
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such  that  bt  is  a column  of  §tt  but  not  of  Btt>  and  such  that  the 
coefficient  of  t>t  in  the  representation  of  b£  in  (3)  is  nonzero. 

It  follows  that  the  matrix  B^t  arrived  at  by  interchanging  columns 
bfc  and  b£  is  nonsingular.  Since  bt  is  not  in  Btt,  b is  not  in 
B.  Let  B'  be  B with  the  columns  b and  b interchanged.  Then 
§'  is  again  square  block  triangular  with  nonsingular  diagonal  blocks. 

Further,  B'  contains  one  more  column  of  B than  did  B.  Therefore 
the  factorization  BF  is  not  minimal. 

Conversely  if  the  factorization  BF  is  not  minimal  there  is  a 
factorization  B'F'  such  that  B Cl  B'  has  more  columns  than  B Cl  b. 

Thus  for  some  t,  Btt  Cl  has  more  columns  than  Btt  Cl  B . Also, 

since  B^t  is  nonsingular  Btt  Cl  B^t  has  full  column  rank.  Hence 
Btt  Cl  b cannot  be  a basis  for  the  column  space  of  ®tt* 

This  completes  the  proof.  Q 

4.4.  Applications  of  Theorem  4.3  : Theorem  4.3  serves  two  useful 
purposes.  Firstly,  it  provides  us  with  a constructive  method  for 
computing  a minimal  B and  F for  any  B:  For  each  t,  find  a linearly 

Independent  subset,  Mtt»  of  columns  of  B^t  that  is  as  large  as  pos- 

sible. Augment  M with  any  other  convenient  columns,  e.g.,  unit 

columns,  to  form  a square  and  nonsingular  matrix  [M  ' N ],  say. 

cc  • tt 

Then  let  B consist  of  the  columns  that  correspond  to  the  Mtt>  together 

with  the  N , suitably  padded  with  zeros.  B will  thus  have 

[M  I N ] as  its  diagonal  blocks.  The  Theorem  tells  us  that  this  is 

tt  • tt  r 

a minimal  B. 
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Secondly,  from  the  proof  of  the  Theorem,  we  have  a means  of 


constructing  a minimal  factorization  from  any  given  square  block 
triangular  factorization,  B = BF,  as  follows.  Suppose  that  the  jth 
column  of  B,  say  b ^ , is  not  in  B.  For  some  t we  can  write 


0 
?t 
3k 

The  j ‘ column  of  F say  f^,  is  the  representation  of  b in  terms 
of  B.  Since  b^  has  zero  in  the  first  t - 1 periods,  it  follows 
by  square  block  triangularity  of  B,  that  f has  the  form 


and  if  f has  a nonzero  component  corresponding  to  any  column  that 


is  in  B but  not  in  B,  we  can  replace  this  column  by  b.,  and  so 


obtain  a 'better'  factorization.  Repeating  this  must  yield  a minimal 


Remark:  As  a final,  remark  regarding  the  structure  of  F in  relation 


(2),  note  that  G is  simply  that  submatrix  of  F whose  columns  are 


indexed  by  the  columns  of  B that  are  not  in  B,  and  whose  rows  are 


indexed  by  those  columns  of  B that  are  not  in  B 


We  shall  say  that  a column  (row)  of  G is  in  period  t if 


the  corresponding  column  of  B(B)  is  in  period  t 


From  the  above,  it  is  then  clear  that  a factorization  BF  is 


minimal  if  f whenever  a column  of  G is  in  period  t (some  t) , it 


has  zeros  in  all  the  rows  in  period  t 


To  conclude  this  section  we  shall  illustrate  with  the  following 


3 period  example 


4.5.  Example 


To  construct  a minimal  factorization  we  find 


is  already  square  and  invertible,  thus  requiring  no  augmentation. 
However,  we  augment  and  with  unit  vectors  to  yield 


0 

1 

0 

0 

0 


1 0 

2 0 

0 0 

-1  1 

1 0 


Note  that  columns  1,  3,  5 of  B are  columns  1,  2,  4 of  B 
respectively,  and  their  columns  1,  3,  5 of  F are  unit  vectors  with 
ones  in  rows  1,  2,  4. 

Also  columns  3 and  5 of  B are  not  in  B and  columns  2 and  4 
of  B are  not  in  B.  The  elements  in  rows  3 and  5 and  in  columns  2 
and  4 of  F yield 


The  rows  of  G correspond  to  columns  3 and  5 of  B which  in  turn 
are  in  periods  2 and  3 respectively,  while  the  columns  of  G correspond 
to  columns  2 and  4 of  B which  in  turn  are  in  periods  1 and  2 respec- 
tively, as  indicated  below. 


Period  Row  of  F 


Period  1 2 

Col  of  F 2 4 
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For  a factorization  to  be  minimal  it  is  necessary  and  sufficient 
that  =0  if  row  i and  column  j belong  to  the  same  period, 

where  G = (g^).  In  this  case  row  1 and  column  2 both  belong  to 
period  2,  and  g ^ * 0.  Therefore  the  factorization  is  minimal. 

i 1 

I J 


5.  Solving  Linear  Equations 

Two  major  steps  in  the  revised  simplex  method  (Beale  [2])  are 
to  solve  the  equations 


By  = a 


where  y is  the  representation  of  the  entering  column,  a,  and 


B tt  = c 


where  t\  is  the  vector  of  prices  used  to  determine  the  next  entering 
column. 


5.1.  Solving  Against  General  Right  Hand  Sides 

Using  the  factorization  B = BF  to  solve  for  example,  the 
system  By  = a,  we  would,  respectively,  solve  the  systems 


Bz  * a 


Fy  * z 


The  system  Bz  **  a is  easy  to  solve  since  B is  square  block 
triangular.  Partition  z and  a to  conform  with  the  partition  of  B: 
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We  would  solve  respectively 


*11  Z1  = 31 
*22  Z2  “ a2  ” *21  Z1 

*KK  Zk  = aK  ~ *K1  Z1  “ ®K2  Z2  “ ' “ “ *K,K-1  ZK-1 

As  we  shall  show  later  on,  it  will  always  be  possible  to  have 
B be  a submatrix  of  A,  the  original  detached  coefficient  matrix. 

Thus,  finding  the  terms  Bgt  zs>  s > t,  in  the  above  right  hand  side 
would  usually  involve  little  work,  since  the  matrices  Bgt  are  usually 
very  sparse,  and  in  fact  often  zero. 


The  system  Fy  ■ z,  however,  is  more  troublesome.  From  (2) 
we  must  solve 


Qy 


Both  G and  H are  computed  matrices,  and  are  typically  between 

20%  and  40%  dense.  Further,  if  G is  k x k,  then  H is  (m  - k)  x k, 

and  for  m much  larger  than  k,  the  storage  requirements  for  H and 

the  work  involved  in  maintaining  H could  become  prohibitive. 

In  the  following  we  shall  show  however  that  it  is  possible  to 

solve  the  system  By  * a without  knowledge  of  H. 

To  simplify  matters  we  shall  assume  that  we  have  suitably  reordered 

T _ 

the  columns  of  B to  be  BQ  and  the  columns  of  B to  be  BP,  and 
refer  to  these  reordered  matrices  as  B and  B,  so  that 


B - B 


C ,)• 


(For  now  we  are  disregarding  the  original  block  triangular  structure 
we  had  an  B and  B.)  Partitioning  B and  B according  to  the  par- 
tition of  F,  we  may  write 


B - (B 


B2) 


B - (B1 


B2) 
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In  short,  to  solve  By  - a 


we  solve  three  sets  of  equations 


Bz  * a 


1 

y “ 


Bw 


„1  1 

= a - B y 


where  w is 


(13) 


0 

2 

y 

_ - 

T 

There  is  a similar  procedure  for  the  system  B tt  * c,  again 
with  three  sets  of  equations: 


To  verify  this  notice  that  the  first  and  last  relations  in  (14)  yield 


T T — T T 

B A * F B A = F 


- ~ 

" T 

2l 

0 

H 

c 

2 

2 

c 

L c J 

and 


B tt  = F 


-T 
B TT 


“ — 

‘ T „T  2' 

U 



G p + H c 

2 

2 

c 

c 

Therefore  the  second  relation  in  (14)  yields 


so  that 


G p + H 


2 

c 


1 

c 


as  required. 

In  sum,  the  equations  (13)  and  (14)  show  that  in  order  to  solve 
the  systems  (4)  and  (5)  we  need  only  be  able  to  solve  systems  involv- 
ing B and  G.  This  means  that  to  execute  the  simplex  method  the 
storage  and  running  costs  are  now  localized  to  the  factorizations  of 
the  K diagonal  blocks  of  B as  well  as  the  factorization  of  G (that 
is,  apart  from  the  pricing  and  minimum  ratio  operations). 


5.2.  Solving  Against  Structured  Right  Hand  Sides 

Section  4.1  dealt  with  the  case  when  the  right  hand  sides  of 
equations  (4)  and  (5)  are  general  vectors.  However,  most  of  the  time, 
these  right  hand  sides  are  structured,  and  taking  advantage  of  this 
can  cut  down  considerably  on  the  work  involved. 
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(1)  In  equation  (4)  the  right  hand  side,  a,  is  almost  always 


an  incoming  column.  Thus  if  the  column  a is  from  period  t,  it 
will  have  the  form 


and  to  solve  the  first  equation  in  (13)  i.e.,  Bz  » a,  we  need  only 
solve  the  reduced  system 


Btt  Zt 


Bt+l,t+l  Zt+1  " at+l  ' Bt+l,t  Zt 


bkk  zk 


*K  ■ bk,i  zi  ' bk,k-i  zk-i 


since  we  know  that 


- 0,  ....  *t-1  - 0 . 
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On  the  average,  we  would  then  solve  [-j]  small  systems,  where 

[•]  denotes  "the  smallest  integer  greater  than."  Note  that  this  would 

not  apply  to  the  third  step  in  (13)  since  the  right  hand  side  there  has 

no  structure  in  general.  Nevertheless  the  total  number  of  small  systems 

K — 

to  be  solved  in  (4)  is  now  on  the  average  K + f-^]  for  B and  1 for 


The  only  time  when  we  are  not  able  to  take  advantage  of  the 
structure  of  a is  when  solving  for  the  current  basic  solution 


Bx  ■ b 


where  b is  the  right  hand  side  of  the  original  L.P.  formulation. 
However,  this  is  only  done  once  every,  say,  50  - 100  iterations. 


(ii)  In  equation  (5)  the  right  hand  side,  c,  is  always  during 
Phase  II,  the  same  unit  vector 


•i- (0 


where  the  objective  row  appears  in  row  Jl  of  the  A matrix.  In 
modern  Implementations  of  the  simplex  method  (Beale  [2]),  this  will 


We  have  assumed  here  that  the  cost  row  (cj_,  C2t  •••»  c^)  in  the 
initial  formulation  has  been  Imbedded  as  one  of  the  rows  of  the  A 
matrix,  as  is  customary. 
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not  be  the  case  during  Phase  I since  c is  set  equal  to  a vector  of 
-l's,  0's  and  l's  depending  on  which  variables  are  currently  outside 
their  bounds. 

For  the  case  c - e^  we  shall  show  that  by  suitably  enlarging 
G,  we  can  solve  for  p in  (14)  without  first  computing  X.  Recall 
that  from  the  derivation  of  equations  (14)  we  had 


However  with  c - e, 


n 

• ^C2J  is  simply  the  column  of  (F^)  ^ . 


Now  since  the  objective  variable  associated  with  the  row  of  A, 

the  cost  row,  is  unrestricted  in  sign,  the  basis  B will  always  con- 
tain its  column,  which  is  also  the  unit  column  e^.  Further  we  can 
always  ensure  that  the  factor  B contains  this  column.  This 
means  that  its  corresponding  column  in  F will  also  be  a unit  vector, 


and  so  not  make  any  contribution  to  G.  By  augmenting  JjjJ  with  this 
unit  vector  we  can  obtain  an  enlarged  G'  with 


c :) 


where  h is  the  corresponding  row  of  H.  Then  the  last  row  of  (G')  * 
will  contain  all  the  nonzeros  of  the  row  of  F“*  and  will  moreover 


- 
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be  precisely  the  p required  in  the  third  equation  of  (14)  . Indeed 

if  we  partition  c =*  e^  into  ( c ^ c2>,  then  c2  *=  0 and  c.  - (0,...,1). 

Thus  in  (14),  A *=  0 during  phase  II,  We  can  therefore  compute 
the  simplex  multipliers,  it,  by  solving  two  systems  of  equations 


and 


(G’)T  p 


-T 
B TT 


(16) 


and  so  achieve  a substantial  savings  over  the  three  systems  in  (14). 

We  conclude  this  section  by  tabulating  the  methods  used  to  solve 
the  systems  (4)  and  (5),  and  also  the  work  involved  in  terms  of  the 
nvmber  of  sma] ler  systems  to  be  solved. 
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6.  Updating  the  Factorization 


Suppose  that  we  begin  an  iteration  of  the  simplex  method  with 
a basis  B and  a minimal  factorization  BF  of  B.  At  the  end  of 
this  iteration  we  obtain  a new  basis  B'  by  replacing  a column  b^ 
of  B with  a column  a of  A.  We  then  need  to  obtain  a minimal 
factorization  B'F'  of  B'  with  as  little  work  as  possible. 

While  we  only  use  the  submatrix  G of  F,  it  will  be  convenient 
here  to  first  determine  how  to  update  f and  then  induce  the  update 
on  G. 

The  two  major  steps  of  the  update  will  be  as  follows: 

(i)  try  to  bring  the  column  a into  B without  disturbing 
square  block  triangularity; 

(ii)  restore  the  factorization  to  minimality  by  seeing  if 
there  are  any  other  columns  of  B,  not  currently  in  B,  that  can  be 
brought  into  B. 

As  we  shall  show,  the  update  on  B can  be  accomplished  by  up- 
dating at  most  two  of  the  factorizations  of  its  diagonal  blocks. 

These  updates  take  the  form  of  an  exchange  of  columns.  The  update  on 

G will  consist  of  a possible  change  in  size,  up  or  down  by  one,  and 

T 

the  addition  of  at  most  two  rank  1 matrices,  each  of  the  form  uv 

T 

where  u is  a column  vector  and  v is  a row  vector. 
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6.1.  Review  of  Updating  an  Inverse  Representation  of  a Matrix  When 
an  Exchange  of  Columns  has  Taken  Place 


Consider  a square  nonsingular  matrix  M,  and  suppose  that  we 
wish  to  modify  M by  replacing  one  of  its  columns  with  some  given 
column  d.  Suppose  further  that  the  columns  we  allow  to  be  replaced 
by  d are  restricted  to  some  submatrix  N of  M. 

Compute  the  representation  y of  d by  solving 

My  *■  d 

If  an  exchange  of  columns  is  at  all  possible  there  will  be  a non-zero 
element  y^  in  y corresponding  to  some  column  m^  of  M that  is 
in  the  restricted  submatrix  N.  Usually,  by  'non-zero'  we  would  mean 
that 

lyJ/'yL  > t°l 


where  TOL  is  a preassigned  tolerance,  typically  of  the  order  10  . 

The  modified  M may  then  be  written  as 


M'  - M + (d  - m1)e^  , 

M'  - H[I  + (y  - ei)e^]  . 
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ombMMB 


\ 


The  inverse  of  the  elementary  matrix  I + (y  - e^e*  is  another  ele- 
mentary matrix  I - 
M may  be  written  as 


1 T 

-Hy  - e.)e  so  that  the  modified  inverse  of 
yi  1 


(M')'1  - [I  - ~ (y  - e)eJ]M  1 . 

yi 

When  the  update  is  actually  implemented  in  this  form,  it  is 
called  'updating  in  product  form'  (see  e.g.,  Dantzig  [6]).  However, 
there  are  today  several  other  methods  available  to  accomplish  this 
update.  These  depend  on  the  representation  of  M ^ as  well  as  the 
degree  to  which  one  wishes  to  preserve  sparsity  and/or  numerical  stabi- 
lity. See  for  example  Forrest  and  Tomlin  [9],  and  Saunders  [23]. 

6.2.  Updating  B 

Let  the  entering  column  a replace  the  i ^ column  b^  of  B. 
Let  a be  in  period  t i.e.,  a has  the  form 


r 

i 0 


a 

t 


*K 
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We  wish  to  determine  whether  a can  replace  a column  of  Bt  where 


5t4 


'tt 


Kt 


From  the  calculation  of  the  representation  of  a in  terms  of  B 
(see  equations  (13))  we  already  have  zt,  the  representation  of  at 


in  terms  of  B 


tt 


Thus  at  can  replace  any  column  of  Btt  having  a 


nonzero  coefficient  in  zfc.  At  this  point  the  following  are  possible. 


6.2.1. 


(a)  The  column  leaving  B is  in  Bt>  and  has  a nonzero 


coefficient  in  z . Column  a then  replaces  b^  in  B by  an  exchange 
of  their  tth  period  components  in  Btt-  (b)  Notice  that  we  have 
changed  §tt  H Btt>  the  basis  for  Btt>  and  it  is  now  possible  that 
®tt  ^ Btt  is  n°  ^on8er  a basis  for  B^.  Thus  to  restore  the  factori- 
zation to  minimality  we  scan  the  columns  of  B^t  that  are  not  in  B^t 
for  a candidate  to  replace  some  column  of  B^t  that  is  not  in  B^. 

This  can  be  done  as  Indicated  in  Section  4.4. 


If  6.2.1  is  not  the  case  then: 
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6.2.2. 

Assuming  that  the  column  a is  not  already  in  B,  we  try  to 
introduce  it  into  Bt  in  place  of  a column  that  is  not  in  B.  If 
there  is  no  suitable  nonzero  in  zfc  to  permit  this,  the  column  a 
cannot  enter  B. 


6.2.3. 

The  column  leaving  B may  be  in  B.  Suppose  it  is  in  period 

s (it  is  possible  that  s = t) . We  then  scan  the  columns  of  B in 

period  s that  are  not  in  B for  a candidate  to  replace  b^^  in 

B . This  can  be  done  by  generating  the  corresponding  row  of  B ^ B 
& ss  s s 

and  then  picking  the  column  corresponding  to  a suitable  nonzero  in 

til  — 

this  row.  Thus  if  b.  is  the  k column  of  B , we  solve 

1 s 


-T 

B w - e. 
ss  k 


T T 

and  the  required  row  is  then  v * w B . If  we  succeed  in  finding 

s s 

T 

a suitable  nonzero  in  v corresponding  to  column  bj,  say  of  B, 
we  then  represent  b^  in  terms  of  B by  solving  B r * bj . 

From  6.2.1,  6.2.2  and  6.2.3  we  see  that  the  update  on  B thus 
takes  the  form  of  an  exchange  of  columns  in  either  B or  Bgg  or 
both,  where  t and  s respectively  are  the  periods  to  which  the  enter- 
ing and  leaving  columns  belong. 


39 


6.3.  Updating  G 


In  updating  G we  must  consider  the  effect  on  F of  the 
column  interchanges  in  B and  in  B.  It  is  best  to  treat  the  cases 
in  6.2  jointly  with  several  subcases.  Combinations  of  the  following 
are  possible: 


I:  a is  column  b^  in  B, 

•k  — 

I : a is  not  a column  of  B, 


II: 

bi 

is  column  b,  of 
k 

B, 

* 

II  : 

bi 

is  not  a 

column  of 

B, 

III: 

a 

replaces 

bi  ("  V 

in  B, 

* 

III  : 

a 

replaces 

b£  in  B 

where  b^  is  not  a column  of  B, 

** 

III  : 

a 

does  not 

enter  B, 

IV: 

bi 

(-  bk) 

is  replaced 

in  B by  column  b^  of  B,  where 

is  not  a 

column  of 

B, 

IV*:  column  bi  (-  bfc)  is  not  replaced  in  B, 

V:  b^,  not  in  B,  is  replaced  by  bj , not  in  B. 


* * 

These  can  occur  as  follows:  (I  , II,  III),  (I  , II,  HI,  V), 
(I*,  II,  III*,  IV*),  (I*,  II,  III*,  IV),  (I*.  II,  III**,  IV*),  (I*.  II, 
III**,  IV),  (I*,  II*,  III**),  (I*.  II*.  Ill*),  (I,  II*),  (I,  II,  IV*), 
(I,  II,  IV). 


We  shall  treat  them  individually. 

I , II,  III:  a replaces  the  itb  column  of  B,  b^,  b^  ■ b^.  Note  that 
the  ith  column  of  F is  unit  vector  e^.  Using  the  approach  of  sec- 
tion 6.1  we  may  write 


AO 


■ I I '• 


B'  = B[I  + (y  - e1)e]’] 

and 

B'  * B [I  + (z  - ek)e^] 

where  the  vectors  y and  z are  the  representations  of  a in  terms 
of  B and  B respectively.  Since  B « BF  and  B'  = B'F'  it  follows 
that 


F'  * [I  - j~{z  - ek)ek]  F [1  + (y  - e^e*]  . 

Multiplying  F into  the  third  factor  and  noting  that  Fy  * z and 
Fe^  = ek  We  ®et 


F'  « [I  - -£-(.*  ~ ek)e£l  [F  + (z  - ek)e^]  . 


This  may  be  conveniently  written  as 


IT  T T IT 

r “ F " Tk  z ek  F " ek(eI  - ek  F)  + ^ 2 ei 


We  now  determine  the  update  on  G by  examining  the  correspond- 
ing rows  and  columns  of  the  terms  on  the  right  hand  side  of  this  rela- 
tion. Firstly,  note  that  since  the  entering  column  a has  replaced 
the  same  column  in  both  B and  B,  the  columns  and  rows  used  by  G' 
are  the  same  as  those  used  by  G.  Secondly,  since  row  k and  column 
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T 1 

where  u and  v are  the  corresponding  subvectors  of  z and 

T k 

e^  F,  respectively. 


Since  z is  already  available  to  us  from  an  intermediate  step 


in  the  solution  of  By  = a,  we  have  u immediately  at  hand.  However, 


X th 

e^  F,  the  k row  of  F,  has  no  components  in  G and  needs  to  be 


computed  from  scratch.  To  do  this  we  solve 


-T 

B w “ e. 


for  w,  the  kth  row  of  B \ and  then  set 


T T _1 
v ■ w B 


where  B consists  of  those  columns  of  B that  are  not  in  B. 


* * 

I , II,  III,  V:  Having  performed  the  first  update  in  (1  , II,  III)  we 


have  found  column  bj  of  B'  to  replace  column  b^  of  B'  where 
both  b 


and  b^  are  in  period  t and  are  not  common  to  B'  and 


B'.  Let  r satisfy 


B'  r - b 


j * 
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Then  we  obtain 


F"  - [I  - ~(r  - e^eJlF' 


i.e. , 


F"  * F'  - ^(r  - e£)e*  F* 


Since  the  introduction  of  into  B'  increases  the  number 


of  common  columns, 

the  size  of  G' 

will 

G'  as 

Gil 

g12 

G'  = 

821 

g22 

where 


g21 

822 


is  the  appropriate  subvector  of  r and  [g' 


12 


T 1 

subvector  of  e^  F*  . Then 


nil  _ f I ^ ry  I Dl^ 

G G11  rfc  g21  *12 


nere  g£2  - r 


g^2]  is  the  appropriate 
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* 

I , 


II*,  III 


** 


In  this  case 


B'  = B 


so  that 


i.e. , 


B'  = B[I  + (y  - ei)e^] 
F'  = F[I  + (y  - e1)e][] 
F'  - F + (z  - fi)e^ 


where  f^  is  the  i1"*1  column  of  F.  Since  is  not  a column  of 

B,  f^  is  not  a unit  column  and  the  above  statement  simply  tells  us  to 
replace  column  i of  F with  z.  The  same  applies  to  G,  i.e.. 


G'  = G + uvT 


where  u and  v are  the  appropriate  subvectors  of  z - f^,  and  e^ 
respectively.  Note  that  the  portion  of  f^  corresponding  to  G is 

a column  of  G and  can  thus  be  found  from  the  representation  of  G. 

* * * 

I , II  , III  : Proceeding  as  in  before  we  see  that 

F’  - [I  - i(.  - e^ej]  [F  + (z  - ft)ej]  ’>  . 
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F'  " F “ z/z  “ e£)(e£  F ~ f £i  ei)  + (e£ 

X 


fi>eI 


where  is  the  i 


th 


. th 


is  the  (£,i)  element 


column  of  F and 
of  F.  Since  we  have  added  column  a to  the  common  columns  of  B 


and  B,  the  size  of  G will  be  reduced  by  1.  This  is  born  out  by 

T 

the  term  (e^  - f^e^  in  the  above  relation.  It  says  that  f^  is 

be  replaced  by  unit  vector  e^.  Note  also  that  the  second  term  in 

the  above  right  hand  side  makes  no  contribution  to  the  i^  column  of 

T T 

F*  since  the  i component  of  e^  F - e^  is  zero. 

Partition  G as 


to 


C11 

g21 

G * 

T 

g22 

g21 

where 


*21 


*22 


is  the  subvector  of  corresponding  to  G,  and  [g^ 

T 
£ 


T 1 

the  subvector  of  e.  F corresponding  to  G.  Then 


Note  that  g22  * 


822J  1s 
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C - Gu  + u 

1 T 

where  u is  the  subvector  of  - — z corresponding  to  G^.  g£j 

zk 

can  be  easily  found  from  the  representation  of  G. 

* **  * * * ** 

I , II,  III  , IV  ; This  is  similar  to  the  case  (I  , II  , III  ) 

except  that  f^  is  now  unit  vector  e^.  Thus 

F*  = F + (z  - ek)e]  , 

and  the  size  of  G will  increase  by  1.  Letting  g denote  the  sub- 

X th 

vector  of  z corresponding  to  G,  and  n the  subvector  of  the  k 
row  of  F corresponding  to  G,  we  get 


x * 

To  obtain  h we  proceed  as  in  (I  , II,  III). 

* ** 

I , II,  III  , IV:  Following  6.2.3  suppose  that  r satisfies  B r ■ 
Then  an  analysis  similar  to  that  of  6.3.1  yields 


F' 


ek)e£]  [F  + (z  - ek)e^]  . 
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Then 


We  can  best  view  this  as  a two  step  update.  Let 


F + (z  - ek)ei 


* 1 T * 

F’  = F - — (r  - ek)ek  F 

k 


£ 

F is  simply  F with  its  i column,  ek  replaced  by  z,  and  reflects 

column  a replacing  b^  in  B.  Since  r is  the  representation  of 

- th  * 

bj  in  terms  of  B,  r is  the  j column  of  both  F and  F . Thus 

the  transition  from  F to  F reflects  b^  replacing  b^  (=  bk) 

in  B.  Note  that  the  jtfl  element  of  ek  F is  rk>  so  that  the  j*"*1 

column  of  F'  becomes  e.  . 

k 

Partition  G as  [G^  g]  where  g is  the  corresponding 

subvector  of  r.  Let  u denote  the  subvector  of  z corresponding 
T T 

to  G,  and  let  k denote  the  subvector  of  efc  F corresponding  to 
G^  Then 


I 


G’  = [Gx  u]  - g[hT,  zk]  . 


From  the  update  on  B we  have  both  r and  z so  that  g,  u,  rk»  and  z^ 

T 

are  immediately  available.  However,  h needs  to  be  computed  as  in 


(I  , II,  III), 


I , II,  III  , IV 


Here  we  obtain 


F'  - [I  - f(z  - e^ej]  [F  + (z  - ek)e^] 
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which  simplifies  to 


F'  = F + (e£  ek)e^  - - e£)ej  F . 

T 

The  term  (e£  - ek)ei  reflects  the  fact  that  after  the  update,  column 

SL  of  B is  a common  column  and  column  k is  not,  whereas  before  the 

T 

update,  the  reverse  was  the  case.  Let  h denote  the  row  of  G cor- 
tti  X 

responding  to  the  k row  of  F,  and  let  v denote  the  subvector 

of  the  row  of  F corresponding  to  G.  Then  the  above  statement 

T T 

means  that  we  need  to  exchange  v and  h in  G before  continuing 
with  the  rank  one  modification.  Thus  if  we  partition  G as 


we  get 


where  u is  the  subvector  of  - — z corresponding  to  the  rows  of 

ZSL 
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We  compute  v as  in  (I  , II,  III)  and  obtain  h from  the 

representation  of  G. 

* * 

I , II,  III  , IV:  We  can  think  of  this  case  as  a continuation  of  the 

previous  case,  viz.  perform  the  update  on  G as  described  above,  and 

then  introduce  b^  into  B'  in  place  of  b^  (=  b^) . This  can  be 

* 

accomplished  in  the  same  way  as  the  case  (I  , II,  III,  V)  with  l 
replaced  by  k. 

* * * * 

I,  II,  IV  : This  is  the  same  as  the  case  (I  , II,  III  , IV  ) with  z 

specialized  to  e„,  and  the  resulting  u therefore  being  zero. 

I,  II,  IV:  Proceed  exactly  as  in  (I  , II,  III  , IV). 

* * * * 

I,  II  : We  may  regard  this  as  a special  case  of  (I  , II  , III  ) with 

z replaced  by  and  the  resulting  u being  zero. 

t 

6.4.  Discussion  and  Summary  of  the  Fundamental  Updating  Operations 
In  Sections  6.2  and  6.3  we  enumerated  the  (many)  possible  up- 
dates that  can  occur.  That  these  are  the  only  possibilities  and  that 
minimality  is  indeed  preserved  at  each  iteration  has  been  claimed  but 
perhaps  not  established  with  full  rigor.  However  these  statements 
can  be  easily  verified  with  the  use  of  Theorem  4.3. 

Presenting  the  update  — and  for  that  matter,  implementing  it  — 
has  been  a laborious  task.  However,  it  is  very  much  this  aspect  of  an 
implementation  of  the  simplex  method  that  determines  its  performance 
relative  to  that  of  other  implementations.  The  updating  procedure 
usually  requires  a significant  time  slice,  and  determines  the  stability 
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and  speed  of  the  remaining  iterations.  Reliable  estimates  of  these 
effects  can  only  be  obtained  by  an  in  depth  analysis. 

A glance  at  the  updates  of  the  previous  sections  shows  that 
there  are  several  operations  that  many  of  them  have  in  common.  These 


are: 


(i) 

(ii) 


INVUPD  — updating  the  factorization  of  a B^, 

-T 

FROW  — solving  the  system  of  equations  B w = e^,  and  then 

T 1 th 

computing  w B to  yield  the  structural  part  of  the  k 


(iii) 

(iv) 

(v) 

(vi) 

(vii) 

(viii) 

(ix) 

(x) 

(xi) 

(xii) 


row  of  F, 

FCOL  — solving  the  system  of  equations  B r 


the  representation  r of  some  column 
GROW  — finding  a row  of  G, 

GCOL  — finding  a column  of  G, 

RSWAP  — replacing  a row  of  G, 

CSWAP  — replacing  a column  of  G, 

ADD  — increasing  the  dimension  of  G 
DEL  — decreasing  the  dimension  of  G 


b to  obtain 
in  terms  of  B, 


by  1, 
by  1, 


RANK1  — performing  a general  rank  one  update  on  G,  that  is 
T 

G'  = G + uv  where  u and  v are  arbitrary  vectors, 


B that  is  not  in  B to  replace 


G + uv 

FINDB  — seeks  a column  of 
column  b^  in  B, 

FIND2  — seeks  two  columns,  one  in  B that  is  not  in  B to 
replace  one  in  B that  is  not  in  B. 


We  may  now  use  these  operations  to  depict  the  whole  updating 
process  in  the  accompanying  flow  diagram  in  Figure  1. 

Note  that  the  GCOL  and  GROW  operations  do  not  appear  in  the 
flow  diagram  since  they  are  used  implicitly  in  the  ADD  and  DEL  opera- 
tions . 
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Figure  1 

The  Update  Flow  Diagram 
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From  the  flow  diagram,  we  obtain  the  following  cross  reference 
table  for  the  11  update  combinations  (indicated  by  the  circled  num- 
bers) and  10  operations.  The  entries  in  the  table  are  the  number  of 
times  an  operation  is  used  in  an  update  combination. 


Table  2 

Cross  Reference  of  Update  Operations  and  Update  Combinations 


Estimating  how  much  work  is  done  in  the  individual  operations  can 
only  be  done  for  a specific  implementation.  To  obtain  an  estimate  of 
the  total  work  per  iteration  we  need  to  estimate  the  probability 
distribution  of  the  11  possible  update  combinations.  While  this  is 
dependent  on  the  specific  problem  being  solved,  a reasonable  degree 
of  uniformity  has  been  observed. 


In  the  sections  following  we  shall  address  these  questions  by 
describing  a specific  implementation  and  by  analyzing  the  statistics 
of  a number  of  test  runs. 


I 
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7 . Implementation 

The  foregoing  details  of  the  factorization  were  all  implemented 
in  an  in-core  Fortran  program  LPBLK  on  the  IBM  370/168  at  the  Stanford 
Linear  Accelerator  Center.  Suitable  existing  software  was  adapted 
wherever  possible,  with  the  remainder  of  the  program  being  written  in 
Mortran  [4]^.  Mortran  was  chosen  for  two  reasons:  Firstly,  it  allows 
one  to  do  structured  programming,  a key  ingredient  of  any  developmental 
work  in  this  area.  Secondly,  it  retains  all  the  important  features  of 
Fortran:  fast  run  times,  ability  to  use  the  Watfiv  compiler,  and 
comparability  with  other  implementations,  few  of  which  are  written 
in  any  language  other  than  Fortran. 

Data  is  input  in  MPS/360  format  [18]  including  simple  upper 
and  lower  bounds  on  any  of  the  variables.  The  code  also  has  the 
capability  of  starting  the  simplex  algorithm  from  any  prespecified 
basis,  input  either  in  MPS/360  format  or  in  the  form  of  a bitmap,  which 
is  a vector  of  n elements,  each  indicating  whether  a variable  is 
basic  or  nonbasic  at  its  upper  or  lower  bound.  Such  a capability 
allows  one  to  make  detailed  comparisons  with  other  codes,  especially 
on  large  problems  where  running  the  problem  from  cold  start  to  opti- 
mality is  both  less  instructive  and  of  prohibitively  large  cost. 


Mortran  was  developed  at  the  Stanford  Linear  Accelerator  Center  by 
J.  Cook  and  L.  Shustek.  It  is  a structured  language  together  with 
a Macro  processor  to  translate  the  language  into  Fortran.  Any 
Fortran  statement  is  a Mortran  statement. 
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Several  of  the  important  routines  in  LPBLK  are  adpatations  of 
routines  written  by  John  Tomlin.  Indeed,  LPBLK  is  modeled  on  his 
experimental  code  LPM1  [28],  which  is  an  implementation  of  the  revised 
simplex  algorithm  using  an  LU  decomposition  of  the  basis  with  standard 
product  form  update. 

As  a standard  for  comparison,  we  chose  another  in-core  Fortran 
code,  MINOS  [20]  developed  by  B.  Murtagh  and  M.  Saunders.  The  linear 
programming  part  of  MINOS  is  a fast  and  stable  implementation  of  the 
revised  simplex  method,  taking  no  advantage  of  any  structure  other 
than  sparsity.  It  was  easy  to  implement  in  LPBLK  the  CHUZR  (choose 
pivot  row)  and  PRICE  (choose  entering  column)  routines  of  MINOS  so 
that  any  differences  in  storage  and  run  time  could  be  accounted  for  by 
the  different  methods  of  representing  the  basis  inverse.  For  the 
main  basis  factorization  and  updating  details  of  MINOS,  see  [23].  The 
routines  of  MINOS  are  all  described  in  [24]. 

7.1.  Factorizing  and  Updating  the  Btt 

Since  the  Btt  are  sparse  and  are  updated  by  an  exchange  of 
columns,  it  was  natural  to  treat  them  as  one  would  the  whole  basis  of 
a general  sparse  linear  program.  Accordingly  we  chose  to  do  an  LU 
factorization  of  §tt  and  to  update  it  using  the  method  of  Forrest 
and  Tomlin  [9]. 

For  a given  Btt>  the  LU  factorization  is  computed  as  follows. 
Permute  B to  the  form 
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where  and  are  lower  triangular  and  D is  the  smallest  such 

matrix.  D is  called  the  "bump."  In  order  to  prevent  fill-in  occurring 
in  F,  reorder  the  rows  and  columns  to  obtain 


E U1  F 


where  is  permuted  to  form  an  upper  triangular  matrix.  The 

problem  is  thus  reduced  to  finding  an  LU  decomposition  of  D.  This 
can  be  done  in  a sparse  manner  by  preordering  the  columns  of  D in 
ascending  order  of  "merit,"  i.e.,  compute 


M - l (r  - 1) 

3 V 


where  r^  is  the  number  of  non-zeros  in  row  i,  and  reorder  the  columns 
of  D so  that 


Mj  < M2  < 


Each  Mj  is  an  estimate  of  the  potential  number  of  nonzeros  that  can 
be  created  by  pivoting  somewhere  in  column  j of  D.  Further  details 
of  this  method  may  be  found  in  [25]. 

Forrest  and  Tomlin  updating  is  aimed  at  preserving  the  sparsity 
of  the  LU  factors.  Although  its  stability  properties  are  not  ideal, 
there  is  a convenient  accuracy  test  (see  [26])  that  can  be  used  to 
decide  whether  or  not  reinversion  is  necessary. 


7.2.  Factorizing  and  Updating  G 

From  Section  6,  we  see  that  whichever  factorization  we  choose 
for  G,  it  will  have  to  be  such  that  rank  one  updates  and  changes  in  size 
can  be  easily  and  stably  accomplished.  Unlike  the  Btt»  G is  a com- 
puted matrix  so  that  the  need  for  stable  methods  is  even  greater  for 
G than  for  the  ®tt* 

In  the  light  of  these  considerations,  a QR  factorization  of  G 
seemed  to  be  the  logical  choice.  Since  G turns  out  to  be  typically 
between  30%  and  50%  dense  there  would  be  little  need  for  sparsity 
considerations,  thus  allowing  the  implementation  to  be  relatively 
straight  forward. 

For  a detailed  analysis  of  the  QR  factorization  and  methods  of 
modifying  the  Q and  R factors,  see  [15].  We  shall  provide  only  an 
outline  of  the  main  procedures. 

The  factors  Q and  R are  defined  to  be  those  orthogonal  and 
upper  triangular  matrices  respectively  that  satisfy 


Q G = R 
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If  G is  nonsingular  then  Q and  R are  unique  up  to  some  signs. 
One  method  of  computing  these  factors  is  to  apply  a sequence  of  ele- 
mentary Householder  transformations  to  G to  reduce  it  to  the  matrix 
R.  The  Householder  transformation  is  an  orthogonal  matrix  of  the 


1 - 2 -j  "J 


where  u^  is  a vector  satisfying  u^  u^  =1.  After  r of  these  trans- 


formations have  been  applied  we  obtain 


P2  P1  G 


p - r 


where  G is  p x p and  the  remaining  (p  - r)  x (p  - r)  bottom 
right  hand  matrix  has  yet  to  be  reduced  to  upper-triangular  form. 
Thus  the  vectors  u^  have  p - j + 1 nonzero  components.  After 
p - 1 applications  we  obtain  the  desired  R on  the  right  hand  side 
and  may  form  Q by  evaluating  the  product 


Q - p. 


P P 
2 1 


from  left  to  right. 


is  ill-conditioned  and  that  it  or  the  whole  factorization  should  be 


recomputed  perhaps  with  a different  choice  of  B.  (See  Section  9.3.) 

When  G is  modified  by  a matrix  of  rank  one,  the  factors  Q 
and  R may  be  updated  by  applying  a sequence  of  plane  rotations. 
These  are  orthogonal  matrices  of  the  form 


where  for  a given  vector  x the  numbers  c and 
that  T applied  to  x annihilates  one  of  x^,  x^ 
x^,  k ^ i,j  unaltered. 

T 

Suppose  that  G'  » G + u v and  Q G = R. 


s are  chosen  so 

and  leaves  elements 


Then 


Q G' 


R + u v 


T 


We  now  apply  rotations  in  the  planes  (p  - 1,  p),  (p  - 2,  p),  • • • > 

(1,  p)  respectively  to  annihilate  the  elements  up_^ in 

this  order.  This  yields 


This  operation  is  called  a backward  sweep . Adding  the  elements  of 
the  second  term  on  the  right  to  the  bottom  row  of  R yields 


A forward  sweep  is  now  performed  to  annihilate  the  last  row  of  R 
to  restore  it  to  upper-triangularity . This  is  accomplished  by  means 
of  rotations  in  the  planes  (1,  p) » (2,  p) (p  - lj  p)  respectively. 


n |j 

so  that 

Note  that  the  order  in  which  the  rotations  are  applied  is  crucial. 
Adding  or  deleting  a row  and  column  may  be  done  similarly,  with  certain 
simplifications. 

At  this  point  we  need  to  decide  how  to  handle  Q.  We  can  either 
evaluate  the  product  of  the  rotations  and  Householder  transformations 
and  store  Q as  an  explicit  dense  square  matrix,  or  maintain  Q in 
product  form.  There  are  significant  drawbacks  and  advantages  associated 
with  both  methods,  and  we  shall  analyze  them  in  detail  in  Section  8. 

7.3.  Recomputing  the  Factorization  B = BF 

The  ability  of  a simplex  code  to  recompute  the  representation 
of  the  basis  in  an  efficient,  sparse  and  stable  manner  is  often  what 
distinguishes  it  from  codes  employing  different  basis  representations. 

In  a typical  run,  one  would  refactorize  the  basis  say  every  50-100 
iterations  in  order  to  reduce  the  number  of  nonzeros  currently  in  the 
representation  and  also  to  retain  numerical  accuracy.  When  starting 
the  run  from  a prespecified  basis,  the  initial  step  is  to  refactorize 
this  basis  before  proceeding. 

In  Section  4 we  sketched  a method  of  finding  a minimal  factori- 
zation for  any  given  B.  This  involved  finding  a maximal  linearly 
independent  subset  of  columns  of  each  B and  then  augmenting  them 
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appropriately  with  unit  columns  to  obtain  the  corresponding  Btt*  It 
turns  out  that  the  method  of  computing  the  LU  factors  of  a sparse  non- 
singular matrix  outlined  in  Section  7 . 1 can  be  easily  adapted  for 
this  task. 

We  begin  by  permuting  B to  the  form 


as  before,  the  only  difference  being  that  D is  now  rectangular.  P 
may  have  more  or  fewer  rows  than  columns,  and  may  be  of  any  rank. 
Proceeding  exactly  as  before  we  compute  the  "merit"  counts  of  each 
column  of  D and  rearrange  them  so  as  to  be  in  ascending  order  of 
merit. 

We  now  simply  pivot  in  this  order,  skipping  a column  when  there 
is  no  nonzero  of  some  suitable  magnitude  available  as  a pivotal  element. 
After  one  such  pass,  a second  pass  can  be  made  through  the  set  of 
rejected  columns  in  case  some  of  them  are  now  suitable  for  pivoting. 
Finally,  we  Insert  unit  columns  in  those  rows  in  which  no  pivoting 
has  yet  taken  place. 

The  "merit"  counts  used  in  this  manner  are  still  reasonably 
good  estimates  of  the  potential  number  of  rionzeros  that  can  be  created 
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by  pivoting  in  a column.  This  is  because  like  the  matrices  Btt» 

D is  usually  not  very  far  from  being  square,  and  is  close  to  being  of 
full  rank. 

It  is  usually  difficult  to  answer  the  question,  what  is  the 
rank  of  a matrix,  or  more  precisely  in  our  context,  what  is  our 
threshold  for  determining  whether  or  not  the  next  column  is  linearly 
dependent  on  the  current  set.  Observe  that  in  our  case,  stopping  a 
few  columns  short  of  obtaining  a maximal  independent  set  is  of  little 
consequence  since  it  only  implies  that  we  shall  be  working  with  a 
factorization  that  is  slightly  off  from  being  minimal.  Thus  in  the 
interests  of  stability  it  would  pay  us  to  use  stricter  tolerances 
than  in  the  case  when  we  know  what  the  rank  of  a matrix  should  be. 

Once  we  have  found  our  factor  B we  represent,  in  terms  of  B, 
each  column  of  B that  was  rejected  in  the  above  process.  Thus  if 
p columns  were  rejected,  we  solve  p sets  of  equations  of  the  type 


By  = b± 


(some  i) 


( 


and  form  G by  extracting  from  each  y those  components  corresponding 
to  the  augmented  unit  vectors. 


7.4.  Looking  at  tne  Fundamental  Operations 

We  are  now  in  a position  to  reconsider  the  factorization  algorithm 
with  respect  to  the  fundamental  operations  outlined  in  Section  6.4. 

From  the  above  implementation,  it  can  be  seen  that  each  of  these  operations 
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can  be  broken  down  into  combinations  of  the  following  basic  operations: 

TRANB  — a transformation  involving  one  of  the  Btt, 

TRANQ  — a transformation  involving  Q, 

TRANR  — a transformation  involving  R, 

SWEEP  — applying  a sequence  of  rotations  to  annihilate  a vector 
of  nonzeros. 

By  considering  the  specific  details  of  the  implementation 
we  arrive  at  the  following  table. 

Table  3 


Cross  Reference  of  Operations  in  the  Update 


TRANQ 


INVUPD 

FROW 

FCOL 

RSWAP 

CSWAP 

ADD 

DEL 

RANK1 

FINDB 

FIND2 


TRANB 

TRANQ 

(Q  EXPLICIT) 

.5 

cfi 

1.5 

1 

.5 

1 

1 

P 

K 

i — , i - -I 
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Each  entry  is  the  average  number  of  times  that  the  operation 


heading  its  column  is  used  by  the  operation  heading  its  row.  For 
example,  the  RSWAP  operation  uses  1.5  TRANQ  operations  when  Q is 
stored  explicitly. 

Note  that  the  TRANR  operation  is  not  used  during  the  update.  It 
is  used  only  twice  per  iteration,  in  the  solution  of  equations  (4)  and 

(5).  Also  from  Table  1 we  see  that  the  TRANB  operation  is  in  addi- 

K K 

tion  used  3K  + [y]  times  during  Phase  I and  2K  + [y]  times  during 

Phase  II. 


r 


\ 


8 . Estimating  the  Work  per  Simplex  Iteration  Accounted  for  by  the 
Factorization 

Obtaining  a good  estimate  of  how  much  work  is  involved  in  any 
implementation  is  difficult  if  not  impossible.  Factors  such  as  the 
machine  and  compiler  in  question,  array  accesses,  the  characteristics 
of  the  problem  being  run  and  the  like,  all  have  an  important  bearing 
on  the  ultimate  runtime.  However,  a criterion  that  is  often  used  by 
numerical  analysts  is  the  total  number  of  multiplications.  While  it 
is  far  from  adequate,  it  is  usually  instructive  to  analyze  an  algorithm 
from  this  point  of  view. 

We  shall  now  do  this  by  estimating  how  many  multiplications  per 
simplex  iteration  are  accounted  for  by  each  of  the  two  factors  B and 
G.  Closely  related  to  estimating  this  quantity  for  G is  the  question 
of  G in  product  form  versus  G stored  explicitly. 

8.1.  Q Stored  in  Product  form  Versus  Q Stored  Explicitly 

A priori  we  see  that  storing  Q in  product  form  saves  us 
evaluating  the  product  of  the  householder  matrices  as  well  as  applying 
each  rotation  explicitly  to  Q.  However  transformations  involving  Q 
become  more  expensive,  and  in  cases  where  a row  or  column  of  Q is 
desired  it  now  has  to  be  computed  by  transforming  an  appropriate  unit 
vector.  Further,  as  the  factors  accumulate,  the  need  may  arise  for  a 
more  frequent  refactorization  of  G. 
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In  the  following  we  shall  assume  that  the  size  of  G remains 
roughly  constant  at  p * p and  that  G and  its  QR  factorization  are 
recomputed  every  q iterations.  We  shall  also  require  the  following 
statistics 

(i)  a : the  average  number  of  sweeps  per  iteration, 

(ii)  p : the  average  proportion  of  rotations  per  sweep, 

(iii)  y : the  average  number  of  transformations  involving  Q,  i.e., 

TRANQ  operations, 

(iv)  n : the  average  number  of  nonzeros  in  the  representation  of 


Note  that  in  working  with  n,  we  are  assuming  that  the  B^'s 
are  all  comparable  in  size  and  further  that  the  frequency  for  recomput- 
ing the  factorization  B = BF  has  been  fixed  beforehand,  so  as  to 
keep  r)  constant. 

Let  us  begin  by  estimating  the  number  of  multiplications  required 
to  generate  the  p columns  of  F of  which  G is  a submatrix.  Solv- 
ing a system  of  the  type 


Solving  systems  of  the  type 


\t  yt ' rhs 


will  require  at  most  ri  multiplications.  For  sparse  right  hand  sides 
one  would  usually  require  significantly  fewer  multiplications.  How- 
ever the  right  hand  sides  become  progressively  more  dense  as  t 

£ 

increases.  On  the  other  hand  l-j]  underestimates  the  number  of  small 
systems  to  be  solved,  so  that  we  may  reasonably  hope  that  the  two 
estimates  off-set  each  other  and  conclude  that,  on  the  average,  the 
number  of  multiplications  to  generate  G is 


r IS.  , 

p n [y]  • 


8.1.1.  Q Explicit 

It  is  straightforward  to  show  that  recomputing  the  QR  factori- 
zation and  forming  the  product  of  the  Householder  transformations  both 
3 2 

require  2/3  p + 0(p  ) miltiplications.  Thus  if  we  recompute  and 
refactorize  G every  q iterations,  we  obtain  a contribution  of 


"7  { P n [f]  + T p3  + 0(p")} 


multiplications . 
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Applying  the  rotation 


( 


to  a vector  of  2 components  requires  4 multiplications.  Since  Q is 

square  and  R is  triangular  we  find  that  a full  sweep  of  pp  rota- 

2 2 
tions  requires  2p  p + 0(p)  multiplications  for  R and  4pp  + 0(p) 

multiplications  for  Q,  yielding  a contribution  of 


6p  p 


2 


+ 0(p) 


multiplications  per  sweep. 

2 2 

Also,  TRANQ  and  TRANR  operations  require  p and  p /2  + 0(p) 
multplications  respectively. 

Noting  that  the  number  of  transformations  per  iteration  is  y 
for  Q and  2 for  R (one  for  the  simplex  multipliers  and  one  for  the 
representation  of  the  incoming  column)  and  that  we  are  assuming  a 
sweeps  per  iteration,  we  obtain  a total  estimate  of 

^ (pn  if]  + y p3^  + (6  a p + y + Dp2  + o(p) 

multiplications  per  iteration,  assuming  that  q ” 0(p).  This  is  the 
total  number  of  multiplications  per  iteration  that  can  be  attributed 
to  G when  Q is  stored  explicitly. 
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8.1.2.  Q In  Product  Form 


3 

Here,  recomputing  Q and  R costs  us  2/3  p + 0(p)  multi- 
plication per  iteration  since  we  do  not  compute  the  product  of  the 
Householder  transformations. 

Each  transformation  involving  Q requires  the  application  to 

a vector  of  p - 1 Householder  transformations  of  diminishing  size, 

2 

a total  of  p + 0(p)  multiplications.  Further,  if  there  are  r rota- 
tions currently  in  the  representation  of  Q,  these  will  require  4r 
multiplications.  Since  each  sweep  produces  pp  rotations,  there 
will  be  appj  rotations  after  j iterations.  The  average  of  this 
quantity  over  q iterations  is  1/2  appq  + 0(p)  rotations.  There- 
fore each  transformation  involving  Q requires,  on  the  average, 

2appq  + p2  + 0(p) 

multiplications. 

As  in  8.1.1  observe  that  a sweep  applied  to  R requires 

2 

2p  + 0(p)  multiplications,  and  that  transformations  involving  R 
2 

require  p /2  + 0(p)  multiplications. 

Doing  a sweeps,  y transformations  with  Q and  2 with  R 
we  finally  obtain  a total  estimate  of 

P T)  [§J  + § P3}  + 2apypq  + (y  + 2a  + l)p2  + 0(p) 
multiplications. 
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8.2. 


The  Contribution  of  B 


At  the  end  of  Section  7.4  we  noted  that  we  require  3K  + 

TRANB  operations  during  Phase  I and  2K  + [-^J  TRANB  operations  during 
Phase  II.  In  addition  we  need  to  consider  the  number  of  TRANB  operations 
used  during  the  update. 

Let  £ denote  the  number  of  FROW  or  FCOL  operations,  and  p 
the  number  of  TRANB  operations  in  the  remaining  update  (INVUPD,  FINDB, 
FIND2) . Then  the  total  number  of  TRANB  operations  is 

3K  + (5  + 1)  [f]  + P 


during  Phase  I and 


2K  + + 1)  [f]  + P 


during  Phase  II. 

Multiplying  these  by  n yields  the  required  number  of  multipli- 
cations. 

8.3.  Summary 

We  summarize  in  the  following  table. 


Table  4 


Average  Total  Number  of  Multiplications  Per  Iteration 


B 

PHASE  I 

{3K  + (£  + 1)  [y]  + p}n 

PHASE  II 

{ 2k  + (C  + l)  [f]  + y)n 

G 

EXPLICIT 

PRODUCT 

uMmaiiiiiiM 

These  estimates  are  difficult  to  interpret  without  knowing 
what  the  values  of  the  constants  are.  In  Section  9 we  shall  use  the 
test  data  to  gain  more  insight  into  this  analysis. 


9.  Computational  Results  and  a Further  Analysis 

As  our  test  problems  we  used  5 dynamic  models.  Their  statistics 
are  summarized  below. 


Table  5 

Problem  Statistics 


Rows 

Structural 

Columns 

626 

1376 

204 

202 

491 

1169 

398 

2750 

301 

480 

% 

Density 


Periods  Non-zeros 


6026 

551 

4029 

11334 

3372 


PIL0T8  is  an  8 period  SIGMA  version  of  the  PILOT  energy  model 
currently  under  study  at  Stanford  University.  It  has  a staricase 
structure  together  with  a few  nonzeros  sprinkled  in  the  lower  block 
triangle.  For  a description  of  the  model  see  [8]. 

The  other  models  are  staircase  models  supplied  by  James  Ho. 

They  are  all  derived  from  real  models  in  economic  planning  (SC205), 
agricultural  production  scheduling  (SCRS8) , engineering  design  (SCSD8) 
and  dynamic  traffic  control  (SCTAP1).  For  more  detailed  descriptions 
of  these  models  see  [13]. 
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The  purposes  of  presenting  the  data  are  twofold:  (i)  to  obtain 
a basic  understanding  of  the  behavior  of  this  factorization  method,  and 
(ii)  to  establish  its  potential  as  a means  of  solving  large  dynamic  time 
period  models.  In  the  former  case  we  recorded  observations  relevant  to 
the  analyses  of  the  preceding  sections,  as  well  as  making  runs  with 
varying  degrees  of  partitioning,  e.g.,  combining  every  3 periods  to  form 
new  periods  each  3 times  as  large  but  1/3  as  many  in  number.  In  addition 
is  recorded  a run  of  PIL0T8.S,  the  PIL0T8  model  rescaled  using  the  geo- 
metric mean.'''  The  scaling  was  done  by  Mohammed  Aganagic. 

The  runs  themselves  were  made  in  several  different  ways.  For 
PIL0T8  and  PIL0T8.S  we  started  from  an  advanced  basis  and  stopped  the  run 
after  2 minutes  CPU  time.  This  yielded  approximately  300  iterations. 

The  same  advanced  basis  was  used  for  all  these  runs.  For  SCSD8,  we 
started  from  an  advanced  basis  and  allowed  the  runs  to  terminate  at  opti- 
mality, again  requiring  about  500  iterations.  The  remaining  models  were  all 
started  from  an  identity  basis  and  allowed  to  terminate  at  optimality. 

The  iteration  counts  for  these  were  230  for  SC205,  870  for  SCRS8,  and 
400  for  SCTAP1. 


1 


I 


I 1/  2 

ror  each  row  compute  (maxla  |min|a  I)  “ where  the  min  is  taken  over 

i 3 j 3 

a^j  j*  0.  Divide  the  row  by  this  quantity.  Proceed  similarly  for  the 
columns.  See  [27]. 


75 


Analysis  of  the  Estimated  Number  of  Multiplications  per  Iteration 


We  began  by  assessing  the  probability  distribution  of  the  eleven 
update  combinations  discussed  in  Section  6.4.  This  was  done  by  record- 
ing the  update  at  each  iteration  and  then  computing  the  relative  fre- 
quencies. These  appear  in  Table  11  at  the  end  of  this  section. 

Following  this  we  used  Table  11  together  with  the  cross  reference 
Table  2 to  obtain  the  average  frequency  per  iteration  with  which  each  of 
the  operations  INVUPD,  FROW,  ...  is  used.  These  figures  appear  in 
Table  12.  The  values  for  y,  a,  and  y in  Table  13. were  computed 
from  Table  12  together  with  the  cross  reference  Table  3.  The  remain- 
ing constants  in  Table  13  were  taken  directly  from  the  individual  runs. 

A glance  at  Table  13  shows  that  the  constants  a,  y,  and  to  some 
extent  y and  p,  are  remarkably  independent  of  the  model,  the  number  of 
periods  (K)  and  the  size  of  G (p) . Indeed  it  seems  reasonable  to 
conclude  that  for  any  problem  similar  in  nature  to  the  ones  tested, 
we  may  take 

5 * i 

y * 2 

a = 2.5 

Y ~ 3.2  Q explicit 

3.7  Q product 

P s .6 

If  we  do  this  and  approximate  [K/2]  by  K/2,  we  obtain  from  Table 
4 the  following  estimates  of  the  number  of  multiplications  involved. 
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Product 


Explicit 


Table  6 

Average  Total  Number  of  Multiplications  per  Iteration 


Phase  I 

(4  K + 2)n 

Phase  II 

(3  K + 2)n 

Explicit 

q-  {I  ? K n + 1 

p3}  + 13.2  p2 

Product 

“ {\  p K n + | 

p3}  + 11.1  p q + 9.7  p2 

The  first  question  to  consider  is  what  is  the  best  value  for 
q,  the  frequency  with  which  we  recompute  G and  its  factorization. 
Plotting  the  number  of  multiplications  for  G as  a function  of  q 
yields  the  following  graph. 

Figure  2 


mults 
for  G 


With  Q stored  explicitly,  it  clearly  never  pays  to  recompute 
G when  the  number  of  multiplications  is  the  only  consideration.  How- 
ever, for  stability  reasons,  it  turns  out  to  be  essential  to  recompute 
G fairly  frequently.  A good  rule  of  thumb  seems  to  be  q s 30. 

With  Q stored  in  product  form,  the  theoretical  minimizer, 
qm  (say),  can  easily  be  shown  to  be 


qm  = .23 J~K  n + P 


The  following  table  gives  the  number  of  multiplications  for 
each  of  the  test  runs.  These  were  found  by  substituting  into  Table 
4 the  observed  values  of  the  constants  in  Table  11. 


\ 

I 


Table  7 

Total  Multiplications/Iteration 


B (1000 ' s) 

G (1000 ’s) 

Total 

Name 

K 

Phase 

qm 

Explicit 

Product 

Phase  II 

I 

II 

q - 30 

q = qm 

Explicit 

PIL0T8 

8 

13.9 

10.5 

23 

81.3 

77.6 

91.8 

4 

20.5 

16.0 

17 

39.7 

40.3 

55.7 

PIL0T8 . S 

8 

14.6 

11.1 

23 

67.8 

69.2 

78.9 

4 

22.7 

18.1 

17 

31.4 

33.8 

49.5 

SC205 

19 

1.2 

.9 

3 

.3 

.6 

1.2 

SCRS8 

4 

8.8 

6.8 

12 

.3 

.8 

7.1 

39 

4.1 

3.1 

12 

32.4 

34.7 

35.5 

20 

6.8 

5.1 

12 

14.9 

17.7 

20.0 

SCSD8 

13 

8.5 

6.5 

11 

7.7 

10.7 

14.2 

8 

10.1 

7.7 

12 

4.4 

6.5 

12.1 

5 

12.5 

9.8 

14 

1.0 

2.1 

10.8 

SCTAP1 

10 

1.9 

1.5 

6 

.4 

.9 

1.9 

5 

2.9 

2.3 

7 

.4 

.8 

2.7 

78 


From  this  table,  we  see  that  with  the  exception  of  the  PIL0T8 
rim  with  K = 8,  the  explicit  form  for  Q is  always  superior.  How- 
ever it  should  be  noted  that  there  will  be  models  for  which  the  product 
form  is  much  superior.  These  will  be  ones  having  high  values  of  p, 
i.e.,  G very  large,  and  high  values  of  p,  the  i.e.,  G very  dense. 
Also  in  cases  where  the  storage  of  Q explicitly  requires  much  core, 
the  product  form  too  will  have  an  advantage. 

As  far  as  the  degree  of  partitioning  is  concerned,  notice  that 
refining  the  partition  has  the  effect  of  increasing  the  size  of  G 
(necessarily)  and  decreasing  the  sizes  of  the  diagonal  blocks  of  B. 
Thus  we  would  expect  some  intermediate  partition  to  be  best.  However, 
we  see  that  on  PIL0T8,  PILOT8.5  and  SCSD8  it  is  far  better  to  have  a 
small  number  of  periods  than  a large  number.  The  main  reason  for  this 
is  that  a disproportionate  amount  of  computation  goes  into  maintaining 
G,  so  that  an  increase  in  the  size  of  G has  a far  greater  effect 
than  a corresponding  increase  in  the  number  of  rows  of  a period.  This 
makes  it  clear  that  any  improvement  in  the  methods  of  representing  G 
will  be  very  worthwhile. 

Finally,  let  us  remark  that  the  number  of  multiplications  has 
- i n^d  out  to  be  a reliable  indicator  of  the  work  involved.  In  Figure 
' -w , . plotted  the  observed  CPU  time  per  iteration  (see  Table  8) 

...I  number  of  multiplications  estimated  in  Table  7. 


Figure  3 

Graph  of  Time/Iteration  Versus  Multiplications/Iteration 


With  the  exception  of  the  4 leftmost  data  points  (corresponding 
to  the  small  problems  SC205,  SCRS8,  and  SCTAP1)  we  see  that  this  is  a 
very  good  straight  line  fit. 
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9.2.  Different  Strategies  During  Update 


So  far  we  have  considered  the  update  and  the  subsequent  analysis 
only  from  the  point  of  view  of  maintaining  the  minimality  of  G at 
every  step.  Maintaining  the  minimality  of  G is  expensive  and  should 
clearly  only  be  done  if  the  effect  of  not  doing  so  is  more  expensive. 

From  the  update  flow  diagram  in  Figure  1 we  see  that  the  points 
at  which  we  do  perform  work  to  maintain  minimality  are  the  FIND2  and 
FINDB  operations.  Let  us  consider  each  of  these  in  turn. 

Recall  that  the  FIND2  operation  scans  the  columns  of  B (in 
some  fixed  period)  that  are  not  in  B for  a candidate  to  replace  a 
column  of  B (in  the  same  period)  that  is  not  in  B.  If  there  are 
l columns  of  B in  this  period  that  are  not  in  B then  failure  to 
find  a suitable  candidate  costs  us  precisely  H TRANB  operations. 

If  we  succeed  in  finding  a candidate,  the  work  involved  is  on  the 
average  i/2  TRANB  operations  and  one  each  of  the  FCOL,  INVUPD,  DEL 
and  RANK1  operations.  What  we  gain  is  a G with  its  size  reduced 
by  one. 

What  will  be  of  great  importance  here  is  the  number  of  itera- 
tions remaining  before  we  recompute  the  whole  factorization  and,  in 
so  doing,  restore  G to  minimality.  Denote  this  number  by  r.  In 
order  to  make  the  following  analysis  possible  we  shall  need  to  assume 
that  if  we  do  not  reduce  the  size  of  G by  one  — when  it  is  possible 
the  remaining  r iterations  will  have  a G that  is  always  one  larger 
than  it  would  have  been.  In  addition,  we  shall  require  an  estimate 
of  the  probability  of  a FIND2  operation  being  successful.  Denote  it 
by  0.  A very  crude  estimate  of  0 may  be  made  by  setting 


Pr{update  combination  5} 

P^.  {update  combination  A}  + P^iupdate  combination  5} 


where  the  probabilities  on  the  right  hand  side  are  given  in  Table  11. 
Refer  to  Figure  1 and  Table  2 for  the  meaning  of  "update  combination 


Now  from  the  cross  reference  Table  3 we  see  that  an  FCOL  opera- 
tion requires  [K/2]  TRANB  operations,  that  an  INVUPD  operation  re- 
quires half  a TRANB  operation,  and  that  the  DEL  and  RANK1  operations 
require  1.5  TRANQ's  and  2.5  SWEEP’S.  Using  arguments  as  in  Section 
8.1.1  this  amounts  to 


uf]  + f)n  + (1.5  + I5p)p2 


multiplications.  Using  a value  of  p = .6  and  approximating  [K/2] 
by  K/2  this  reduces  to 


|(K  + l)n  + 10.5  p2 


Let  W(p)  denote  the  work  per  iteration  (in  multiplications) 
when  G is  p x p.  Then  assuming  that  the  size  of  G stays  roughly 
constant,  the  expected  net  profit  for  a FIND2  operation  is 


0{r[W(p)  - W(p  - 1)1-  ±(K  + l + l)n  - 10.5  p2}  - (1  - 9Hn  • 


From  Table  11  we  obtain  for  the  explicit  form,  with  q = 30 


W(  p)  - W(  p - 1)  ~ ■—  K n + p2  + 26  p . 

Taking  as  an  example  the  PIL0T8  model  with  K = 4,  p = 46,  n 51  1117, 
l - 12  (obtained  by  guessing  at  i = p/K)  and  9 = .05  (estimated  as 
suggested  above)  we  obtain  the  expected  profit  to  be 

14320  - 78r  . 

This  is  positive  for  r > 184. 

We  therefore  conclude  that  with  p = 46  and  fewer  than  184 
iterations  to  go  before  a complete  refactorization,  it  does  not  pay 
us  to  use  the  FIND2  operation.  However  on  a problem  like  SCSD8  with 
K =*  39,  p=  44,  n “ 26,  £ = 2 (say),  and  0 = .02  it  turns  out  that 
the  critical  value  of  r is  17.  Here  it  pays  us  to  u^e  the  F1ND2 
operation  most  of  the  time. 

With  the  F1NDB  operation,  we  may  proceed  similarly.  Here  we 
seek  a column  of  B that  is  not  in  B to  replace  a given  column  of 
B.  Going  back  to  the  update  flow  diagram  in  Figure  1 we  see  that  the 
FINDB  operation  occurs  in  three  positions.  For  update  combinations  2, 
3 and  8,  9 we  see  that  we  either  fail  to  find  a suitable  column  and 
do  no  extra  work,  or  we  succeed  and  proceed  with  one  each  of  the  FCOL, 
INVUPD,  DEL  and  RANK1  operations,  exactly  as  in  the  FIND2  case.  The 


: 
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only  difference  between  this  and  the  FIND2  case  is  that  whether  we 
find  a column  or  not,  there  is  initially  only  one  TRANB  operation. 

Thus  the  expected  net  profit  becomes 

0{r[W(p)  - W(p  - 1)]  - -|(K  + l)n  - 10.5  p2}-n  . 

For  the  PIL0T8  model  with  the  same  parameters  as  above,  except  for  a 
new  estimate  of  0 - .25  (obtained  as  in  the  F1ND2  case)  we  obtain 
the  critical  value  of  r to  be  17.  For  the  same  SCSD8  model  as  before, 
we  get  0 = .38  and  a critical  value  of  r = 40.  The  remaining  FINDB 
operation  in  the  flow  diagram  may  be  handled  likewise. 

Using  ideas  of  this  type  for  fine  tuning  will  be  valuable, 
especially  on  models  where  many  runs  are  made,  and  reasonable  estimates 
of  the  various  parameters  are  available. 

9.3.  A Comparison  with  MINOS 

We  ran  each  of  the  problems  on  MINOS,  adjusting  tolerances  and 
refactorization  frequencies  where  necessary  to  obtain  run  times  that 
were  as  fast  as  possible.  The  runs  on  LPBLK  were  all  made  with  Q 
stored  explicitly,  and  maintaining  the  minimality  of  G at  each  step. 
The  following  table  resulted.  (See  over.) 

ETA  is  the  number  of  nonzeros  in  the  representation  of 
each  Btt  for  LPBLK  and  the  whole  basis  for  MINOS.  The  times  are 
given  in  CPU  milliseconds.  All  figures  were  computed  as  averages  over 
the  whole  run. 
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LPBLK 


MINOS 


Name 


PILOT8 


PILOT8 . S 


SC205 


SCRS8 


SCSD8 


SCTAP1 


K 


19 


39 

20 

13 

8 

5 


10 

5 


ETA(n) 


423 

1117 


433 

1168 


15 


502 


47 

130 


Size  of 
G (p) 


71 

46 


66 

42 


26 

— 

mm 

82 

n 

Wmm 

131 

163 

295 

17 

163 

545 

7 

158 

6 

5 


Time/ 

Iteration 


302 

240 


288 

235 


30 


92 


55 

56 


ETA 


14582 


15815 


1021 


3335 


4358 


2204 


Spikes 


152 


130 


23 


38 


56 


28 


Time/ 

Iteration 


203 


188 


25 


65 


130 


40 


From  the  CPU  times  in  the  above  table,  we  see  that  MINOS  is 
approximately  22%  faster  than  LPBLK  on  PILOT8,  PIL0T8.S,  SC205  and 
SCSD8.  It  is  much  faster  on  SCRS8  and  SCTAP1. 

Scaling  seems  to  have  had  the  same  effect  on  both  codes.  In 
the  case  of  MINOS  the  scaled  columns  would  make  the  pivot  elements  of 


triangle  columns  more  acceptable  (in  a relative  test)  and  so  force 
fewer  triangle  columns  to  become  spikes.  It  would  appear  that  a simi- 
lar argument  could  be  made  for  LPBLK:  if  more  pivot  elements  become 

, I acceptable  then  the  FINDB  and  FIND2  operations  would  be  more  succes- 

sful  in  keeping  the  size  of  G to  a minimum. 
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On  the  stability  side,  the  only  models  to  present  problems  were 


PIL0T8  and  PIL0T8.S.  In  the  runs  recorded  here,  LPBLK  experienced 
no  trouble;  however  there  were  some  advanced  bases  for  PIL0T8  that 
LPBLK  failed  to  refactorize  while  MINOS  succeeded.  As  mentioned  in 
Section  7.2  the  ratio 


is  a lower  bound  on  the  condition  number  of  G.  On  PIL0T8  and  PIL0T8.S 

g 

this  ratio  was  typically  of  the  order  10  , sometimes  going  as  high 

10  3 

as  10  . With  the  other  models,  this  ratio  was  at  most  10  . The 

PILOT  model  has  always  been  badly  conditioned  and  a code  like  MPSIII 
[19]  often  needs  to  force  unit  columns  into  the  basis  and  lose  feasi- 
bility in  order  to  maintain  stability.  MINOS  is  perhaps  the  only 
code  that  has  not  suffered  unduly  on  PILOT. 

One  apparent  cause  of  PILOT'S  illconditioning  is  the  presence 
of  dense  12  x 12  Leontief  matrices  imbedded  in  each  period.  The 
variables  corresponding  to  these  Leontief  systems  were  all  free^  (no 
upper  or  lower  bounds)  and  hence  always  in  the  basis.  In  the  case 
of  LPBLK,  by  the  manner  in  which  the  Btt  are  defined  during  a complete 
refactorization  (see  Section  7.3),  the  dense  Leontief  columns  would 
be  considered  last  for  introduction  into  B and  many  of  them  would 


This  was  done  since  it  was  possible  to  demonstrate  in  advance  that 
in  any  optimal  solution,  these  activities  would  always  be  at  a posi- 
tive level. 


not  'make  it.'  This  would  result  in  several  columns  of  G being 
transforms  of  these  Leontief  columns,  a fact  that  seemed  intuitively 
wrong.  Accordingly  we  forced  all  the  Leontief  columns  into  B by 
setting  their  "merit"  counts  artificially  at  zero.  Having  then  chosen 
each  B in  this  way,  we  would  refactorize  the  ®tt's  immediately 
to  achieve  sparsity,  after  lifting  the  restrictions  on  the  "merit" 
counts.  The  result  of  doing  this  was  a much  improved  factorization: 
the  lower  bound  on  the  condition  number  of  G dropped  by  as  much  as 
two  orders  of  magnitude. 

9.4.  Storage  Considerations 

Table  9 below  contains  statistics  on  the  average  number  of  non- 
zeros required  to  represent  the  basis  factorization.  The  first  two 
columns  give  the  average  number  of  nonzeros  (for  both  LPBLK  and  MINOS) 
taken  over  the  whole  run.^  The  third  column  gives  the  average  number 
of  nonzeros  recorded  immediately  after  a MINOS  refactorization. 

This  table  shows  that  the  basis  storage  requirements  of  LPBLK 
are  between  23%  and  60%  of  those  of  MINOS,  a substantial  savings. 

The  method  of  Bisschop  and  Meeraus  [3]  leaves  the  initial  factor 
ization  of  the  basis  untouched  and  instead  updates  a certain  augmented 
matrix  of  size  k x k where  k is  the  number  of  columns  in  which  the 
initial  and  current  bases  differ.  A lower  bound  on  their  storage 
requirements  is  therefore  given  by  the  number  of  nonzeros  in  the 


^These  were  taken  directly  from  Table  8,  with  the  entries  in  the  LPBLK 
column  being  set  equal  to  K n + 3/2  p2  for  the  explicit  form  of  Q. 
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Name 


10946 

7642 


14582 


15815 


Initial  representation.  Assuming  that  their  initial  factorization 
is  an  LU  as  performed  by  MINOS  we  may  use  the  third  column  of  Table 
9 as  the  lower  bound.  Comparing  this  with  the  column  for  LPBLK,  we 
observe,  at  the  appropriate  degree  of  partitioning  that  LPBLK  is  far 
superior  on  all  except  SCRS8  and  SCSD8  where  it  is  no  worse  that  the 
lower  bound . 

Table  9 

Storage  Requirements 


Another  possibility  is  to  use  the  factorization  of  LPBLK  as  the 
initial  factorization  in  the  Bisschop-Meeraus  scheme.  However,  while 
this  will  more  or  less  equalize  the  storage  requirements  it  seems 
very  unattractive  from  a stability  point  of  view. 


Average  Over  The  Run  Average  at 

refactorization 

LPBLK  MINOS  (MINOS) 


PIL0T8 


11874 


PIL0T8.S 


12088 


SC205 


SCRS8 


SCSD8 


SCTAP1 


9.5.  Bumps  and  Spikes 


It  is  interesting  to  compare  the  bump  and  spike  structure  of 
a basis  with  its  natural  block  triangular  structure.  Perhaps  surpris- 
ingly the  two  appear  to  have  little  in  common. 

Table  10  below  gives  the  bump  and  spike  structure  of  a typical 

PIL0T8  basis  (626  * 626).  This  was  found  using  Hellerman  and  Rarick's 
4 

p algorithm  as  implemented  in  MINOS. 

Table  10 

Bump  and  Spike  Structure  of  a PIL0T8  Basis 


BUMP// 

//COLUMNS 

//SPIKES 

1 

2 

1 

2 

3 

2 

3 

2 

1 

4 

2 

1 

5 

2 

1 

6 

2 

1 

7 

2 

1 

8 

2 

1 

r\ 

1 

K9- 

2 

1 

13 

8 

381 

100 

Ka 

3 

2 

Notice  that  there  are  several  very  small  bumps  and  a single  very  large 
one.  The  block  triangular  structure  on  the  other  hand  contains  8 
blocks  each  having  approximately  80  rows.  Thus  in  such  a case,  a 
method  like  that  of  McBride  [17]  is  unlikely  to  be  effective. 
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Probability  Distribution  of  the  11  Update  Combinations 


Average  Number  of  Times  an  Operation  is  Used  During  Update 


Table  13 


Name 

K 

€ 

'■  1 ■ 

n 

■ 

Y 

HI 

u 

H 

Ot 

Explicit 

Product 

■ 

PIL0T8 

8 

.75 

1.86 

a 

2.23 

1 

3.81 

.55 

4 

.88 

2.60 

□ 

2.30 

3.89 

.74 

PIL0T8.S 

8 

.83 

2.42 

■ 

66 

2.46 

3.83 

.47 

4 

.91 

3.65 

M|||| 

42 

2.35 

EuS 

3.79 

.66 

SC205 

19 

1.20 

.93 

15 

4 

2.70 

3.35 

3.88 

.95 

SC..S8 

4 

1.08 

1.30 

502 

4 

2.49 

3.21 

3.59 

.38 

39 

BS9 

1.23 

26 

2.97 

3.39 

4.11 

.52 

20 

1.54 

82 

m 

2.75 

3.30 

3.85 

.51 

SCSD8 

13 

1.07 

1.67 

155 

19 

2.66 

3.24 

3.73 

.61 

8 

1.06 

2.03 

295 

17 

2.49 

3.20 

3.58 

.50 

5 

1.06 

1.79 

545 

7 

2.30 

3.11 

3.30 

.55 

SCTAP1 

1.10 

tm 

47 

6 

2.41 

3.53 

.42 

5 

1.11 

130 

5 

2.42 

BH 

3.52 

.51 

Key  to  Table  13 
K : number  of  periods 

£ : average  number  of  times  per  iteration  that  an  FROW  or  FCOL  opera- 
tion is  used 

y : average  number  of  TRANB  operations  per  iteration  other  than  those 
in  FROW  and  FCOL 

0 : average  number  of  nonzeros  in  the  representation  of  a typical 

Btt. 

p : average  size  of  G 

a : average  number  of  SWEEP  operations  per  iteration 
Y : average  number  of  TRANQ  operations  per  iteration 
p : average  proportion  of  rotations  formed  per  sweep. 
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10.  Conclusion 

As  matters  stand,  the  computational  results  of  the  previous 
section  show  that  LPBLK  saves  considerably  on  storage  requirements. 
However  from  the  run  times  it  is  apparent  that  LPBLK  comes  close  to, 
but  is  not  yet  competitive  with  the  best  available  software.  It  seems 
that  problems  on  which  it  will  excel  will  be  ones  like  PILOT 8 with 
more  periods  and  perhaps  a greater  proportion  of  non-zeros  in  the 
lower  block  triangle. 

However  the  preceding  analysis  clearly  shows  that  cheaper  ways 

have  to  be  found  to  deal  with  G.  Any  method  that  can  halve  the  work 

involved  with  G will  make  a 30%  difference  in  the  run  time,  so  making 

it  very  competitive  on  the  test  problems  presented.  One  possibility 

T 

currently  under  study  is  to  use  the  factorization  Q G = R but  then 

discarding  Q altogether  and  instead  maintaining  G explicitly.  To 

solve  the  equations  Gx  = g we  would  first  solve  two  triangular 
T T 

systems  R Ry  = g,  and  then  set  x = G y.  Since  G is  typically 

between  30%  and  50%  dense  it  may  well  be  possible  to  store  it  as  a 

T 

sparse  matrix,  so  that  forming  G y will  be  inexpensive.  Thus  solv- 
ing equations  this  way  will  involve  no  more  work  than  when  we  keep  Q. 
Where  we  save  considerably  is  in  the  update,  since  now  we  need  only 
apply  the  rotations  to  R. 

This  approach  was  first  used  by  Gill  and  Murray  [10],  in  the 

context  of  the  simplex  method  for  dense  systems.  Paige  [21]  has 

shown  that  the  above  method  of  solving  the  equations  is  stable,  involv- 

2 

lng  only  K(G)  (instead  of  K (G))  in  the  error  bound.  However, 
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to  solve  for  the  price  vector,  we  need  to  solve  equations  of  the  form 
T 

G x = g (refer  to  Section  5.1).  This  requires  the  solution  of  the 
T 2 

equations  R Rx  ■ Gg  where  now  the  term  K (G)  does  enter  into  the 

error  bound.  Since  the  price  vector  only  Indicates  which  column  to 

bring  into  the  basis,  and  is  not  used  in  any  further  computation,  less 

stability  here  can  be  tolerated.  Further,  in  Phase  II  we  can  simplify 

T 

the  procedure  since  we  always  solve  G x * ef  where  er  is  a fixed 
unit  vector.  Multiplying  this  relation  on  the  left  by  Q yields 
Rx  = Qer*  Therefore  by  maintaining  the  rC**  column  of  Q explicitly 
from  one  iteration  to  the  next,  we  can  solve  these  equations  by  means 
of  a single  triangular  system.  This  is  both  stable  and  inexpensive. 

In  addition  to  the  choice  of  representation  for  G,  the  fine 
tuning  aspects  of  Section  9.2  have  yet  to  be  properly  Implemented  and 
tested.  Some  experimental  runs  made  along  these  lines  have  been  very 
encouraging . 

More  work  needs  also  to  be  done  on  the  stability  side  in  general. 
There  is  a definite  need  for  more  care  in  deciding  which  columns  go 
into  B and  which  are  transformed  to  become  columns  of  G. 


for  his  careful  reading  of  the  manuscript  and  helpful  suggestions. 
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Time— staged  and  multi-staged  linear  programs  usually  have 
a structure  that  is  block  triangular.  Basic  solutions  to  such  prob- 
lems typically  have  the  property  that  similar  type  activities  persist 
in  the  basis  over  several  consecutive  time-periods.  When  this 
occurs  the  basis  is  close  to  being  square  block  triangular.  In 
1955  Dantzig  suggested  a way  of  factorizing  the  basis  to  take  advan- 
tage of  this  property.  This  paper  discusses  persistence  in  stair- 
case models  and  then  presents  a considerable  refinement  of  the  above 
factorization  algorithm. 

The  method  has  been  implemented  in  an  experimental  code, 
with  use  being  made  of  L V and  QR  factorization  and  updating  tech- 
niques for  the  solution  of  small  sub-systems  of  equations . An  in- 
depth  analysis  is  made  of  the  work  involved,  and  computational  exper- 
ience on  several  dynamic  models  is  reported. 
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